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Prefaces 


Digital image processing is an important research area. The techniques developed in this area so far 
require to be summarized in an appropriate way. In this book, the fundamental theories of these techniques 
will be introduced. Particularly, their applications in image denoising, restoration, and segmentation will 
be introduced. The entire book consists of four chapters, which will be subsequently introduced. 


In Chapter 1, basic concepts in digital image processing are described. Chapter 2 will see the details of 
image transform and spatial filtering schemes. Chapter 3 introduces the filtering strategies in the 
frequency domain, followed by a review of image restoration approaches described in Chapter 4. 
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1. Introduction 


1.1 Digital image processing 


Digital image processing is the technology of applying a number of computer algorithms to process digital 
images. The outcomes of this process can be either images or a set of representative characteristics or 
properties of the original images. The applications of digital image processing have been commonly found 
in robotics/intelligent systems, medical imaging, remote sensing, photography and forensics. For example, 
Figure 1 illustrates a real cardiovascular ultrasonic image and its enhanced result using a Wiener filter 
that reduces the speckle noise for a higher signal-to-noise ratio. 


(a) (b) 


Figure 1 Illustration of image enhancement by applying a 
Wiener filter to a cardiovascular ultrasonic image: (a) original, 
(b) enhanced image. 


Digital image processing directly deals with an image, which is composed of many image points. These 
image points, also namely pixels, are of spatial coordinates that indicate the position of the points in the 
image, and intensity (or gray level) values. A colorful image accompanies higher dimensional information 
than a gray image, as red, green and blue values are typically used in different combinations to reproduce 
colors of the image in the real world. 


The RGB color model used in the color representation is based on the human perception that has attracted 
intensive studies with a long history. One example area of the RGB decomposition can be found in Figure 
2. The present RGB color model is based on the Young-Helmholtz theory of trichromatic color vision, 
which was developed by Thomas Young and Herman Helmhotz in the early to mid nineteenth century. 
The Young-Helmholtz theory later led to the creation of James Clerk Maxwell’s color triangle theory 
presented in 1860. More details on this topic can be found in [1]. 
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(c) (d) 


Figure 2 An RGB image along with its R, G and B components: 
(a) RGB image, (b) R component, (c) G component and (d) 
B component. 


1.2 Purpose of digital image processing 


The main purpose of digital image processing is to allow human beings to obtain an image of high quality 
or descriptive characteristics of the original image. In addition, unlike the human visual system, which is 
capable of adapting itself to various circumstances, imaging machines or sensors are reluctant to 
automatically capture “meaningful” targets. For example, these sensory systems cannot discriminate 
between a human subject and the background without the implementation of an intelligent algorithm. 


Figure 3 denotes a successful example where a human object is separated from his background using a k- 
means clustering algorithm, which is part of the technologies developed in digital image processing. We use 
this example to justify the importance and necessity of digital image processing. To separate the human 
object from the background, subsequent processes will be employed. These processes can be low-, mid- and 
high-level. Low-level processes are related to those primitive operators such as image enhancement, and 
mid-level processes will get involved in image segmentation, and object classification. Finally, high-level 
processes are intended to find certain objects, which correspond to the pre-requisite targets [2]. The 
following description provides more details about the object segmentation, shown in Figure 3. 
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(a) (b) 


Figure 3 Illustration of human object segmentation from the 
background using a k-means clustering approach [5]. 


Low-level processes are not desirable in this particular example. This is due to the fact that the original 
image has been acquired in a good condition and hence no evidence of image blurring occurs. During the 
mid-level processes, a manual assignment of cluster centers is achieved by a professional. This leads to 
optimal clustering of image intensities by the automatic k-means clustering. If necessary, the extracted 
object will be used to generate human identity, which is one of the high-level processes. Up to now, it is 
clear that without digital image processing one will not be able to generate “meaningful” object in this 
example and beyond. 


1.3 Application areas that use digital image processing 


The applications of digital image processing have been tremendously wide so that it is hard to provide a 
comprehensive cover in this book. While being categorized according to the electromagnetism energy 
spectrum [2], the areas of the application of digital image processing here are summarized in relation to 
the service purpose. This is motivated by the fact that one particular application (e.g. surveillance) may get 
different sensors involved and hence presents confusing information in the categorization. In general, the 
fields that use digital image processing techniques can be divided into photography, remote sensing, 
medical imaging, forensics, transportation and military application but not limited to. 


Photography: This is a process of generating pictures by using chemical or electronic sensor to record 
what is observed. Photography has gained popular interests from public and professionals in particular 
communities. For example, artists record natural or man-made objects using cameras for expressing their 
emotion or feelings. Scientists have used photography to study human movements (Figure 4). 
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Figure 4 Study of human motion (Eakins Thomas, 
1844-1916). 


Remote sensing: This is a technology of employing remote sensors to gather information about the Earth. 
Usually the techniques used to obtain the information depend on electromagnetic radiation, force fields, or 
acoustic energy that can be detected by cameras, radiometers, lasers, radar systems, sonar, seismographs, 
thermal meters, etc. 


Figure 5 illustrates a remote sensing image collected by a NASA satellite from space. 


Figure 5 An example of remote sensing images 
(image courtesy of NASA, USA). 


Medical imaging: This is a technology that can be used to generate images of a human body (or part of 
it). These images are then processed or analyzed by experts, who provide clinical prescription based on 
their observations. Ultrasonic, X-ray, Computerized Tomography (CT) and Magnetic Resonance Imaging 
(MRI) are quite often seen in daily life, though different sensory systems are individually applied. Figure 
6 shows some image examples of these systems. 
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(a) (b) (c) 


Figure 6 Medical image examples: (a) ultrasound, (b) X-ray and (c) MRI. 


Forensics: This is the application of different sciences and technologies in the interests of legal systems. 
The purpose of digital image processing in this field is used to be against criminals or malicious activities. 
For example, suspicious fingerprints are commonly compared to the templates stored in the databases (see 
Figure 7). On the other hand, DNA residuals left behind by the criminals can be corresponding to their 
counterparts saved in the DNA banks. 
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Figure 7 A fingerprint image. 


Transportation: This is a new area that has just been developed in recent years. One of the key 
technological progresses is the design of automatically driven vehicles, where imaging systems play a 
vital role in path planning, obstacle avoidance and servo control. Digital image processing has also found 
its applications in traffic control and transportation planning, etc. 


Military: This area has been overwhelmingly studied recently. Existing applications consist of object 
detection, tracking and three dimensional reconstructions of territory, etc. For example, a human body or 
any subject producing heat can be detected in night time using infrared imaging sensors (Figure 8). This 
technique has been commonly used in the battle fields. Another example is that three dimensional 
recovery of a target is used to find its correspondence to the template stored in the database before this 
target is destroyed by a missile. 


Figure 8 An infrared image of tanks 
(image courtesy of Michigan 
Technological University, USA). 
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1.4 Components of an image processing system 


An image processing system can consist of a light source that illuminates the scene, a sensor system (e.g. 
CCD camera), a frame grabber that can be used to collect images and a computer that stores necessary 
software packages to process the collected images. Some I/O interfaces, the computer screen and output 
devices (e.g. printers) can be included as well. Figure 9 denotes an image processing system. 


Illumination 
Camera 


Output device 
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Frame grabber 


Computer 


Figure 9 Components of an image processing system. 


1.5 Visual perception 


Digital image processing is performed in order that an image does fit the human visual judgments. Before 
our study goes further, the human visual system has to be studied so that the target of image processing 
can be properly defined. In this subsection, we briefly describe the vision principle of the human eye. 


Blood nerve 
wessels 


Figure 10 Illustration of eye anatomy 
(image courtesy of [3]). 
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The reflected light from the object enters the external layer that coats the front of the eye, which is called 
the cornea (see Figure 10). Afterwards, the light passes through some watery fluid namely the aqueous 
humor. Then the light reaches the iris that may contract or dilate so as to limit or increase the amount of 
light that strikes the eye. Passing through the iris, the light now arrives at the pupil, a black dot in the 
centre of the eye, and then reaches the lens. The lens can change its shape in order to obtain focus on 
reflected light from the nearer or further objects. 


1.6 Image acquisition 


Images are generated from the combination of an illuminant source and the reflection of energy from the 
source. In general, images can be two dimensional (2-D) or three-dimensional (3-D), depending on the 
used sensors and methodologies. For example, a set of 2-D cardiovascular images can be piled up to form 
a 3-D image using an automatic correspondence algorithm. These 3-D reconstruction techniques will be 
detailed in later chapters. 


Image acquisition can be categorized to single sensor, sensor strips, and sensor arrays based. For example, a 
photodiode is made of a single sensor (see Figure 11). Computerized tomography uses sensor strips to 

measure the absorption of x-ray that penetrates the human body. Figure 12 shows one of these systems. An 
ordinary camera (e.g. Olympus E-620) is based on rugged arrays that normally consist of millions elements. 
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Figure 11 Illustration of a 
photodiode. 


Figure 12 Illustration of a computerized tomography 
system (image courtesy of [4]). 


1.7 Image sampling and quantization 


1.7.1 Basic concepts in sampling and quantization 


An image consists of an indefinite number of points with continuous coordinates and amplitudes. To 
convert this image to digital form, both the image coordinates and amplitudes must be discretized, where 
the image points will be changed to pixels while the amplitudes use discrete values. Sampling refers to 
digitization of the coordinates, and quantization is the digitization of the amplitude values. 


Figure 13 shows an example of image sampling and quantization, where (b) reveals that the image 
resolution is reduced, whilst (c) indicates that gray levels decrease, compared to (a). 
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Figure 13 Illustration of image sampling and quantization: (a) original, (6) sampling and 
(c) quantization. 


1.7.2 Representing digital images 

An image is normally represented in matrix form, originating from the upper left corner. Also, an image 
consists of a certain number of gray levels. For example, if it has 2™ gray levels, this image is referred to 
“m-bit image”. One of the regular manipulations on the image is zooming or shrinking. Either way may 
bring down the resolutions of the original image due to interpolation or extrapolation of image points. 
One of the commonly noticed image problem is aliasing. This phenomenon appears when the interval of 


image sampling is higher than a half of the distance between two neighboring images points. As a result, 
the Moiré pattern effect will appear and seriously deteriorate the image quality. 


1.8 Basic relationship between pixels 


1.8.1 Neighbors of a pixel 


Assuming that a pixel has the coordinates (x, y), we then have its horizontal and vertical neighbors which 
have coordinates as follows 


(x+1, y), (x-1, y), (x, y+1) and (x, y-1) (1.8.1) 
We can have four diagonal neighbors of the point (x, y) below 

(x+1, y+1), (x+1, y-1), (x-1, y+1) and (x-1, y-1) (1.8.2) 
1.8.2. Adjacency, connectivity, regions, and boundaries 


If an image point falls in the neighborhood of one particular pixel, we then call this image point as the 
adjacency of that pixel. Normally, there are two types of adjacency, namely 4- and 8-adjacency: 


(1) 4-adjacency: Two pixels are 4-adjacent if one of them is in the set with four pixels. 
(2) 8-adjacency: Two pixels are 8-adjacent if one of them is in the set with eight pixels. One 
example is shown below, where red “1s” form the 8-adjacency set. 
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anal 0 0 


Connectivity is relevant but has certain difference from adjacency. Two pixels from a subset G are 
connected if and only if a path linking them also connects all the pixels within G. In addition, G is a region 
of the image as it is a connected subset. A boundary is a group of pixels that have one or more neighbors 
that do not belong to the group. 


1.8.3 Distance measures 

Assuming there are two image points with coordinates (x, y) and (u, v). A distance measure is normally 
conducted for evaluating how close these two pixels are and how they are related. A number of distance 
measurements have been commonly used for this purpose, e.g. Euclidean distance. Examples of them will 


be introduced as follows. 


The Euclidean distance between two 2-D points I(x),y,) and J(x,y2) is defined as: 


f(x, -— =)? +, 2)" (1.8.3) 
The City-block distance between two 2-D points (x;,y1) and (x2, y2 ) can be calculated as follows: 

|x, -x, (+1 -y2 | (1.8.4) 
For the above two 2-D image points, the Chessboard distance is 

max(] x, — x, |, ¥) — V2 |) (1.8.5) 


Summary 


In this chapter, basic concepts and principles of digital image processing have been introduced. In general, 
this chapter started from the introduction of digital image processing, followed by a summary of different 
applications of digital image processing. Afterwards, the fundamental components of an image processing 
systems were presented. In addition, some commonly used techniques have been summarized. 


In spite of their incomplete stories and terse descriptions, these presented contents are representative and 
descriptive. The knowledge underlying with this introduction will be used in later chapters for better 
understanding of advanced technologies developed in this field. 
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Problems 


(1) What is digital image processing? 

(2) Why do we need digital image processing? 

(3) Please summarize the applications of digital image processing in the society. 
(4) How is the human visual perception formed? 

(5) What is the difference between image sampling and quantization? 

(6) What is the difference between adjacency and connectivity? 

(7) How to compute the Euclidean distance between two pixels? 
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2. Intensity transformations and spatial filtering 


2.1 Preliminaries 


In this chapter, we will discuss about some basic image processing techniques that include the intensity 
transforms and spatial filtering. All the techniques introduced here will be performed in the spatial domain. 


2.2 Basic Intensity Transformation Functions 

Photographic Negative: 

This is perhaps the simplest intensity transform. Supposing that we have a grey level image in the 
range[0,1], it is expected to transform the black points (Os) into the white ones (1s), and the white pixels 


(1s) into the black ones (0s). This simple transform can be denoted by (assume that f(x, y) is normalized 
into the range [0,1 ]) 


~ 


f(~%y)=1- f(y) (2.2.1) 


For a 256 gray level image, the transform can be accomplished by 


f(x,y) =1- f(x, y)/256 (2.2.2) 


An example of the Photographic Negative transform is shown in Figure 14. 


(a) (b) 


Figure 14 A panda Image (a), and its photographic negative transformed image (b). 
Gamma Transform 


Gamma transform is also called power-law transform. It is mathematically denoted as follows: 


f(xy=e* f(xy)’ 023) 
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where c and y are two constants. The gamma transform can make pixels look brighter or darker depending 
on the value of y. When f(x, y) is within the range [0,1] and y is larger than one, it makes the image 
darker. When y is smaller than one, it makes the image look brighter. Figure 15 shows the output of the 
Gamma transform against different inputs with the parameters set as 0.5, 1 and 2 (c = 1). From this plot, 
we can see that the curve with y = 2 is below the curve with y = 1. This indicates that the output is smaller 
than the input, which explains why an image in the Gamma transform, when y = 2, will become darker. 
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Figure 15 The plot of the Gamma transform against 
different parameters. Horizontal axis is the input, and 
vertical is the output. 


Figure 16 shows that the Gamma transformed panda image under different parameter values (c = 1). 


(a) (b) 


Figure 16 Gamma transformed images different parameters. (a) Gamma = 2 ; (6) Gamma = 0.5. 
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2.3 Histogram Processing 


Histogram is one simple but very important statistical feature of an image. It has been commonly used in 
image processing. Intensity histogram is a distribution of the grey level values of all the pixels within the 
image. Each bin in the histogram represents the number of pixels whose intensity values fall in that 
particular bin. A 256 grey level histogram is often used, where each grey level corresponds to one bin. 
Using bi to represent the i*” number of bins, the histogram can be represented as follows: 


h(i) = #{@, y), fy) € bi} (2.3.1) 


where # represent the Cardinality of the set. Though it is possible to show the histogram of a color image, in 
this chapter we mainly focus on the gray level histogram. The color histogram will be shown in Chapter 1 
(Digital Image Processing - Part Two). 
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Figure 17 shows an image and its grey level histogram of 64 bins. 


(a) (b) 


Figure 17 A mountain image (a) and its 64 bins grey level histogram (b). 
2.3.1 Histogram Equalization 


It is a fact that the histogram of an intensity image lies within a limited data range. Those images usually 
have black or white foreground and background. Figure 18 shows an image example whose intensity 
distribution is either black or white. From Figure 18 (b), we can see that a very large portion of pixels 
whose intensity rests within the range [0, 50] or [180, 255]. A very small portion of pixel resides in the 
range of [50, 180]. This makes some details of the image hardly visible, e.g. the trees on the mountains in 
the image shown in Figure 18 (a). This problem can be solved by a histogram stretching technique called 
histogram equalization. 


The basic idea of histogram equalization is to find the intensity transform such that the histogram of the 
transformed image is uniform. Of the existing probabilistic theories, there exists such an intensity 


transform. Suppose that we have an image f(x, y), and its histogram h(i), we have the accumulative 
function of h(i) as follows: 


c(i) = f, h(t) dt (2.3.1.1) 


It can be proved that such a transform makes the variable y = c(i) follow a uniform distribution. Thus, 
for a 256 grey level image, the histogram equalization can be performed by the following equation: 


t= a c(f(xy)) (2.3.1.2) 


where n is the total number of pixels in the image. 


Download free ebooks at bookboon.com 


24 


Please click the advert 


Digital Image Processing — Part I 


Intensity transformations and spatial filtering 


(a) 


100 250 400 


Figure 18 (a) The mountain image after histogram equalization, and (b) histogram. 


2.3.2 Histogram Matching 


Histogram Equalization is a special case of histogram matching. The general purpose of histogram 
matching is to find a transform that can adjust an image by converting its histogram to a specified 
histogram a priori. As we have learnt from the above discussion, the desired histogram in histogram 
equalization is a uniform distribution. In this case, there exists a transformation as shown in Eq. (2.3.1.2). 
A straightforward implementation can be achieved. However, in more general case we may want the 
desired histogram to be a specific shape and to adjust the image based on this histogram. This operation is 
called histogram matching. It can be achieved by matching the accumulative distributions of the original 


histogram and desired histogram. Step by 
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1) | Compute the histogram of the original image f (x, y) and its accumulated distribution c(i). 

2) Compute the accumulated distribution of the desired histogram d (i). 

3) | Compute the desired intensity g(x, y) at location (x, y) by matching the d(i) and c(i) using 
minimum distance criteria (nearest n) as follows: 


g(x,y) =i" = argmin,(l\d@ — c(f(x,y))| (2.3.2.1) 


Figure 19 shows some internal results of the above steps on the image shown in Figure 17 (a). 


(c) (d) 


Figure 19 (a) The accumulated histogram of the image in Figure 2.3.1 (a); (b) the 
desired histogram; (c) the desired accumulated histogram; (d) the adjusted image 
based on image matching according to Eq. (2.3.2.1). 
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2.4 Fundamentals of Spatial Filtering 


Filtering is a fundamental operation in image processing. It can be used for image enhancement, noise 
reduction, edge detection, and sharpening. The concept of filtering has been applied in the frequency 
domain, where it rejects some frequency components while accepting others. In the spatial domain, 
filtering is a pixel neighborhood operation. Commonly used spatial filtering techniques include median 
filtering, average filtering, Gaussian filtering, etc. The filtering function sometimes is called filter mask, or 
filter kernel. They can be broadly classified into two different categories: linear filtering and order-statistic 
filters. In the case of linear filtering, the operation can be accomplished by convolution, i.e., the value of 
any given pixel in the output image is represented by the weighted sum of the pixel values of its 
neighborhood (a linear combination) in the input image. For order-statistic filter, the value of a given pixel 
in the output image is represented by some statistics within its support neighborhood in the original image, 
such as the median filter. Those filters are normally non-linear and cannot be easily implemented in the 
frequency domain. However, the common elements of a filter are (1) A neighbourhood (2) an operation on 
the neighbourhood including the pixel itself. Typically the neighbourhood is a rectangular of different 
size, for example 3x3, 5x5 ... and smaller than the image itself. 


2.5 Smoothing Spatial Filters 


2.5.1 Smoothing Linear Filters 


As stated above, linear filters have observed their output values as the linear combination of the inputs. 
The commonly used smoothing filters are averaging and Gaussian filters. The smoothing filters usually 
have the effect of noise reduction (also called low pass filtering, which will be discussed in the next 
chapter). It can be performed using the convolution operation. 


5%, Y) = De Myo Dn nye hen, n) f (x — m,y — n) (2.5.1.1) 


h(m, n) is a filtering mask of size MxN. Each element in this filter mask usually represents the weights 
used in the linear combination. The types of different linear filters correspond to different filter masks. 
The sum of all the elements in h(m, n) will affect the overall intensities of the output image. Therefore, it 
is sensible to normalize h(m,7n) such that the sum is equal to one. In the case of negative values in 
h(m,n), the sum is typically set to zero. For smoothing (low pass filters), the values in the mask must be 
positive; otherwise, this mask will contain some edges (sharpening filters). 


The linear filter can also be performed using correlation as follows: 
SY) = Ln mye Danny h(n n)f (x +m, y +n) (2.5.1.2) 
If h(m,n) is symmetric, i.e. h(m,n) = h(—m, —n), the correlation is equivalent to convolution. 


Average Filtering 


The average filtering is also called mean filtering, where the output pixel value is the mean of its 
neighborhood. Thus, the filtering mask is as follows (3x3 as an example) 
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111 
h=1/9]1 1 1 (2.5.1.2) 
111 


In practice, the size of the mask controls the degree of smoothing and loss of the details. 


Gaussian Filtering 


Gaussian filtering is another important filter. The weights in the filter mask are of a Gaussian function 
(here we assume it is isotropic): 


h(m,n) = = exp(—0.5 «(m2 +n?) /o) (2.5.1.3) 


where a is the variance and controls the degree of smoothing. The larger value a is, the larger degree 
smoothing can be achieved. For example, the 3x3 Gaussian filter with o set as 1 is represented by the 
following equation. 


0.60 0.10 0.60 
0.10 0.16 0.10 
0.60 0.10 0.60 


h(m,n) = (2.5.1.4) 


Figure 20 shows the Gaussian mask of size 5x5, 


Figure 20 The 3D plot of the Gaussian filter with the 
variance set as 1. 


Figure 21 shows the Lenna image as well as the filtered images using mean filtering with different mask sizes 
and Gaussian filtering with different 0. We can see that with the increment of the mask size of the mean filter, 
the level of smoothing becomes stronger. This also holds for the Gaussian filtering with different o. 
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(9) 


Figure 21 (a) The Lenna image; (b) (c) (d) filtered images using mean filtering with mask size 3, 
7, 11; (e) (f) (g) filtered images using Gaussian filtering with different variances at 1, 5, 9. 


Download free ebooks at bookboon.com 


29 


Digital Image Processing — Part I Intensity transformations and spatial filtering 


Please click the advert 


2.5.2 Order-Statistic Filters 


The order-statistical filters are usually non linear filters, which are hardly represented by convolution. 
Commonly used filters include median filter. There are some other filters as well such as max/min filter. The 
median filter simply replaces the value of the pixel with the median value within its neighborhood. Similarly, 
the max/min filter replaces the value of the pixel with the maximum/minimum value within its neighborhood. 


Figure 22 shows the filtered image using a median filter of different mask sizes. 
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2.6 Sharpening filters 


The objective of sharpening is to draw the attention to the fine details of an image. This is also related to 
the situation where an image that has been blurred and now needs to be de-blurred. In contrast to the 
process of image smoothing that normally uses pixel averaging techniques, sharpening can be conducted 
using spatial differentiation. The image differentiation actually enhances edges and other discontinuities 
and depresses the areas of slowly changing gray-level values. 


The derivative of a function is determined via differencing. A basic first-order derivative of a one- 
dimensional function f(x) is 


F _ px41)- f(a) 2.6.1) 
ox 


Similarly, a second-order derivative can be defined as follows: 


arf 


Ox? ~ 


F(xt)+ f(x-D-2f (x) (2.6.2) 


In this section, we mainly talk about the use of two-dimensional and second-order derivatives for image 
sharpening. It is inevitable to mention the concept of isotropic filters, whose response is independent of 
the discontinuity direction in the image. 


One of the simple isotropic filters is the Laplacian, which is defined as follows with a function f(x,y): 


Intensity transformations and spatial filtering 


2 2 
Vf= fy f (2.6.3) 
Ox" Oy 
Where 
a’ f 
a hE Ay) 2 y) (2.6.4) 
aa 
By =f(x,y+)+ f(x, y-l)-2f (x,y) (2.6.5) 
Then the two-dimensional Laplacian is obtained summing the two components: 
Vf =[f@+Ly)+ f(x-Ly)+ fy +D+ fy-DI-4 f(y) (2.6.6) 
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To implement Eq. (2.6.6) we illustrate four filter masks used for the Laplacian in Figure 23. 


0 1 0 
1 -4 1 
0 1 0 
1 1 1 
1 -8 1 
1 1 1 
0 af 0 
4 4 | 
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1 “ “4 
“4 8 aj 
“1 4 4 


Figure 23 Filter masks used in the Laplacian enhancement. 
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Such a Laplacian enhancement can be described as 


fee ne VFN) oe 
F(x, y)+ Vo fy) 
To simplify the Laplacian we have 
F(x, y)= f(x y)-[fathLy)+ f(e-Ly)t+ fa ytD+ f%y—D] (2.6.8) 


One of the examples using the Laplacian is shown in Figure 24. 


1000; 


Figure 24 Image sharpening and histograms: (a) original image and its histogram 
(c), (b) shapened image and its histogram (d). 
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For a function f(x, y), the gradient of fis defined as 


f 
Vf = s (2.6.9) 
Oy 
In the implementation, this equation can be changed to 
Vi = (2, +22, +2,)—(@, +223 +2,)|+| 4, +22, 4+2,) — (2, +22, +2,) | (2.6.10) 


Where (z;,..29) are the elements of the 3x3 mask. This mask can be one of the following patterns: 


Z4 Zo 23 
Z4 Z5 Z6 
Z7 Zg Z9 
1 e) z 
0 0 0 
1 2 1 
-1 0 1 
-2 0 2 
A 0 1 


Figure 25 Image mask that is used to generate gradients. 


2.7 Combining image enhancement methods 


In real applications, it is hard to know what kind of noise has been added to an image. Therefore, it is 
difficult to find a unique filter that can appropriately enhance this noisy image. However, it is possible if 
several de-blurring methods can be combined in a framework in order to pursue a maximum denoising 
outcome. In the section, we explain one of the combinatorial techniques using an X-ray example [6]. 


Figure 26 (a) illustrates a male chest’s X-ray image. The purpose of the process is to highlight the middle 
cross section of the image using the combination of sharpening and edge detection. Figure 26 (b) shows 
the result of applying the median filter, (c) is the outcome of using Sobel edge detection, and finally (d) 
demonstrates the combination of the Laplacian and Sobel process. Figure 26 (d) shows the edge details in 
the graph, which highlights the structure of the central cross section of the image. 
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e— 


(c) (a) 


Figure 26 Illustration of using the combination of different segmentation algorithm: (a) original, (b) 
median filtering, (c) Sobel edge detection and (d) combination of Laplacian and Sobel process. 


Summary 


In this chapter, some basic image intensity transformation functions have been provided. This consists of 
negative and Gamma transforms. Afterwards, histogram based processing techniques were introduced. 
Especially, histogram equalisation and matching techniques were summarised due to their importance in 


image processing applications. 
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To further enhance the image quality we then introduced a number of spatial filters. Particularly, we 
discussed about averaging filtering, Gaussian filtering, and statistic filters. Due to the needs of enhancing 
image contrasts, we then focused on the technical details of sharpening filters, which have been commonly 
used in the image processing applications. As an extension, we discussed about the applications of the 
combination of image enhancement methods. One of the examples is to integrate a Laplacian filter with 
the Sobel edge detection for better image details. 


References 


[6]  http://www.radiologyinfo.org/en/photocat/gallery3.cfm?pid=1&Image=chest- 
xray.jpg&pg=chestrad, accessed on 15 October, 2009. 
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Problems 


(8) Please perform the image negative and Gamma transform for the cameraman image as follows: 


(9) Please use the Matlab tools to generate codes to add Gaussian noise (mean=0, variance = 2) to 
the above cameraman image, followed by averaging and Gaussian filtering respectively. 
(10) Can you find a proper sharpening filter for the following motion blurred image? 


(11) Why do we need to combine different filtering approaches for image enhancement? 
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3. Filtering in the Frequency Domain 


3.1 Background 


The Fourier transform is a fundamental mathematical tool in image processing, especially in image 
filtering. The name of ‘Fourier transform’ or ‘Fourier Series’ can be traced back to the date of year 1822 
when a French mathematician published his famous work - the Analytic Theory of Heat. In this work, he 
had shown that any periodic function can be expressed as the sum of sines and cosines of different 
frequencies weighted by different values. This ‘sum’ interpretation is usually termed as ‘Fourier series’ (in 
modern mathematical theory this is in fact a special case of orthogonal decomposition of a function within 
one of the functional space e.g. Hilbert space, where each axis is represented by a function, e.g., the sine 
or cosines functions.). For non-periodic function, this type of operations is called ‘Fourier Transform’. For 
simplicity, in the following sections we will use the general term ‘Fourier transform’ to represent such a 
transform in both of the periodic and non-periodic functions. 


The merit of the Fourier transform is that a band limited function or signal could be reconstructed without 
loss of any information. As we will see, this in fact allows us to work completely in the frequency domain 
and then return to the original domain. Working in the frequency domain is more intuitive and can enable 
the design of image filters easier. 


Although the initial idea of the Fourier transform was applied to the heat diffusion, this idea has been 
quickly spread into other industrial fields and academic disciplines. With the advent of computer and 
discovery of the fast Fourier transform, the real-time practical processing of the Fourier transform 


becomes possible. 


In the following sections, we will show step by step the basic concept of Fourier transform from 1-D 
continuous function to the 2-D case. 


3.2 Preliminaries 


3.2.1 Impulse function and its properties 


Impulse function is one of the key concepts of sampling, either in spatial domain or in the frequency 
domain. A typical impulse function is defined as follows: 


co ifx=0 
= 2.1 
5) {0 otherwise Re) 
And its integral is subject to the following constrains: 
fo 6@)dx =1 (3.2.2) 
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A very important feature of the impulse function is that for any continuous function f(x), the integral of 
the product of the function f(x) with the impulse function is simply the value of the f(x) at the location 
where the impulse happens. For example, if the location of the impulse is at T, the following mathematical 
expression can help us to get a sample from the function f(x) at T. 


Jo FCS (x —T)dx = f(T) (3.2.3) 


The above equation implies that the operation of getting a value from the function f(x) at T can be 
performed by the integral of the product of this function with the impulse function at location T. This is 
the fundamental of the sampling theory. 


3.3 Sampling and the Fourier Transform of Sampled Functions 


3.3.1 Sampling 


In the real world, a scene or an object could be mathematically represented by a continuous function. However, 
in today’s digital world, to show or store an image of the scene in a computer the scene has to be converted into 
a discrete function. This digitalization process in fact can be called ‘sampling’. In fact, sampling is the 
process/operation of converting a continuous function into a discrete function. One role of this operation is that 
it makes the image representation tractable in a computer memory or storage and displayable on screen as well 
as visible to humans. In this process, a key question is what kind of sampling frequencies should be adopted in 
the sampling process (or which point will be selected to describe the discrete version of the function)? A 
sample refers to a value or set of values at a point in time and/or space. Mathematically, the sampling principle 
can be described as follows (beginning from simple periodic functions): 


What do you see? 


A skateboarder? 


We see an All Options Junior Trader leaving 
the office after a successful day. 
He has the freedom to make his own decisions. 
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Consider a continuous 1-D periodic function, y(t) = y(t + nT) with n = 0,1,2,3, ..., where t represents 
the time in seconds. Now we are going to sample this function by measuring the value of continuous 
function every T seconds. T here is called sampling interval. Thus the sampled function (discrete version) 
can be represented by 


yln] = y(nT),n = 0,1,2,3, .... (3.3.1) 


The sampling interval is an indicator of the frequency that we use to sample the function. Here we 
introduce a sampling frequency or rate that measures the number of samples in unit (usually in one 
second), f = 1/T. The sampling rate is conventionally measured in hertz or in samples per second. 
Intuitively, an ideal sampling rate allows us to reconstruct the original function successfully. 


The output of sampling is a numerical sequence of the continuous function, usually represented by a 
matrix or vector of numbers. Image can be viewed as a 2D sampling version of the real world. A 
conventional representation of the image is a matrix of M number of rows and N number of columns. The 
spatial coordinates denotes the relative location of the pixels within the matrix (element index) while the 
entry of the matrix represents the value of the grey level intensity or color intensity at that pixel, usually 
represented by I (i, j). The smallest coordinates are conventionally set as (1,1) in the popular MATLAB 
program. Note the notation (i, j) in a matrix does not correspond to the actual physical coordinates. For 
example, a function f(x); x = 0,0.5,1 is represented by a matrix [f (0), f (0.5), f(1)]. The entry value at 
(1,2) is f (0.5) with actual physical coordinates x at 0.5. 


3.3.2 The Fourier Transform of Sampled Functions 


The convolution theory tells us that the Fourier transform of the convolution of two functions is the 
product of the transforms of the two functions. Thus by applying the Fourier transform on the both side of 
Equation 3.2-3, we can get: 


F(f(nT)) = F(F(O@d(0) 
- eae F(u-=) G22) 


This shows that Fourier transform of the sample function is a sampled version of the Fourier transform 
F(u) of continuous function f(t). Note that the sample pace in the frequency domain is the inverse of the 
sample pace in the spatial domain. 


3.3.3 The Sampling Theorem 


As stated before, the reason why we discuss the sampling and Fourier transforms is to find a suitable sampling 
rate. Insufficient sampling may lose partial information of a continuous function. So, what is the relationship 
between the proper sampling interval and the function frequency? In other words, under what circumstances is 
it possible to reconstruct the original signal successfully (perfect reconstruction)? 
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A partial answer is provided by the Nyquist-Shannon sampling theorem [7], which provides a sufficient 
(but not always necessary) condition under which perfect reconstruction is possible. The sampling theorem 
guarantees that band limited signals (i.e., signals which have a maximum frequency) can be reconstructed 
perfectly from their sampled version, if the sampling rate is more than twice the maximum frequency. 
Reconstruction in this case can be achieved using the Whittaker-Shannon interpolation formula [8]. 


The frequency equivalent to one-half of the sampling rate is therefore a bound on the highest frequency 
that can be unambiguously represented by the sampled signal. This frequency (half the sampling rate) is 
called the Nyquist frequency of the sampling system. Frequencies above the Nyquist frequency can be 
observed in the sampled signal, but their frequency is ambiguous. This ambiguity is called aliasing later 
shown in Figure 27. To handle this problem effectively, most analog signals are filtered with an anti- 
aliasing filter (usually a low-pass filter with cutoff near the Nyquist frequency) before their conversion to 
the sampled discrete representation. We will show this principle using the following examples: 
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Let us consider a simple function f (x) = cos(2mf,x) + cos (27f2x) with f, set as 5 and f2 set as 10. We 
expect that the power spectrum of the FFT of the signal will have peaks at 5 and 10. Figure 27 (a) shows 
the original signal and (b) show the power spectrum. As expected, we can observe two peaks at u = 5 and 
10 respectively in (b). According the Nyquist-Shannon sampling theorem, the minimum sampling rate is 
at least equal to the two times of the highest frequency of the signal. In our example, the trade-off rate is 
20 (2x10). We now sample this signal with a rate larger than 20, say 40 to see what happens. Figure 27 
(c) shows the sampled signal using this rate, and (d) shows the corresponding power spectrum. We can see 
that the signal looks quite similar to the original one with the peaks at frequency 5 and 10. This is 
sometimes called oversampling and leads to a perfect reconstruction of the signal. Another case is to 
sample the signal at a lower rate than 20, for instance 10. Figure 27 (d) and (f) show the sampled signal 
using rate 10 and its power spectrum. From this figure, we can see that trend of the signal is not preserved 
and the power spectrum looks different from the original one in (b). This is called under-sampling or 
frequency aliasing. In practice, most signals are band infinite. That means they usually have a very infinite 
frequency. The effect of aliasing almost exists in every case. However, in the case of high quality imaging 
system such an aliasing is not visually observable by human eyes. 
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Figure 27 Illustration of the sampling theorem: (a) 
original; (b) the power spectrum; (c) the sampled signal at 
sampling rate 40 (oversampling); (d) the power spectrum 
of (c); (e) the sampled signal at sampling rate 10 (under 
sampling); (f) the power spectrum of (e). 


3.4 Discrete Fourier Transform 


3.4.1 One Dimensional DFT 


Let us start with the 1-D discrete Fourier transform (DFT) for a discrete function f (x) forx = 
0,1,2 ...M — 1. Its 1D Fourier transform can be calculated as follows: 


F(u) = Yio f(x) e, (3.4.1) 
w= 2n (=) (3.4.2) 


where u = 0,1,2,...M — 1. wu is usually called frequency variables, similar to the coordinates x in the 
spatial domain. The exponential on the right side of the equation can be expanded into sinuses and cosines 
with variables u determining their frequencies as follows: 
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e J® = cos(w) + jsin(w), (3.4.3) 


The inverse discrete Fourier transform is: 
f(x) = SUI Fae’, (3.4.4) 
3.4.2 Relationship between the Sampling and Frequency Intervals 


Consider a discrete function with the fixed sample interval AT (i.e., the time duration between samples) 
and it contains N samples. Thus the total duration of the function is D = N * AT. Using Au to represent 
the spacing in the Fourier frequency domain, the maximum frequency (i.e., the entire frequency range, 
also the maximum samples that will get during a time unit, usually per second) is just the inverse of the 
sampling interval. 


U = 1/AT (3.4.5) 


Note that we have N data samples and the frequency domain also has the same number of samples. 
1 
AT*N 


The spacing per sample can be calculated as Au = ~ = . Thus, the relationship between AT and Au 


can be stated using the following equation: 


AT = N*Au (3.4.6) 
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3.5 Extension to Functions of Two Variables 


Similar to the 1D Fourier transform, the 2-D discrete Fourier transform can be expressed as follows: 


ux *Y) 


PGi) = Fre 6 deen f (x, ye 2" Carr , (3.5.1) 
Where I(x, y) is the original image of size M x N. u and v are in the same range. 


The inverse Fourier Transform is defined as follows: 


, Ux vy 
I(x,y) = — ying Duco Fu, velar) (3.5.2) 


where F (u, v) is the Fourier transform of the original image. Usually the Fourier transform is complex in 
general. It can be written in the polar form: 


F(u,v) =a+bj = |F(u,v)|efPer (3.5.3) 
And we have the magnitude and the phase angle as follows: 

IFGi, 2) | = (a? +67) (3.5.4) 

o(u,v) = tan), (3.5.5) 


The power spectrum is the square of the magnitude|F (u, v)|. Figure 28 shows an image, its power 
spectrum and phase angle respectively. 


Figure 28 An image of sunflower (a), its power spectrum (b), and its phase component (c). 
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3.6 Some Properties of the 2-D Discrete Fourier Transform 
The 2D Fourier transform has the following properties that are often used. 


1) ‘Translation invariant (magnitude). If a function I(x, y) is translated by an offset (Ax, Ay), its Fourier 

— jam (Ua%4vay a . . 

q 2m MN ). We can see that the translation in the spatial domain 
— jon( tae 242) 

only affects the phase component by e MN 


transform becomes F (u, v Je 


. The magnitude |F (u, v)| is not affected. 
2) Periodicity. The 2D Fourier transform and its inverse are infinitely periodic in the u 

and v directions. 
3) symmetry: 


a. The Fourier transform of a real function is conjugate symmetric: F*(u, v) = F(—u, —v), and the 
magnetite is always symmetric, i.e, |F(u, v)| = |F(—u, —v)|. 
b. The Fourier transform of a real and even function is symmetric: F(u,v) = F(—u, —v). 


4) Linearity: The Fourier transform of the sum of two signals equals to the sum of the Fourier 
transform of those two signals indepedently. g(x,y) + f(x,y) <=> F(u,v) + G(u,v) 

5) Scaling: If the signal is spatially scaled wider or smaller, the corresponding Fourier transform 
will be scaled smaller or wider. The magnitude is also smaller or wider. f (ax, By) <=> 


1 Uv 
ap) G)@): 


The Fourier transform also has other properties. For a complete list of properties, interested readers may 
refer to some handbooks of signal/image processing and sites like Wikipedia and Mathworld. 


3.7 The Basics of Filtering in the Frequency Domain 


As we have stated in Chapter 2, in the spatial domain the filtering operation is sometimes performed by the 
convolution operation with a filtering mask (Note this is the case for linear filter. For nonlinear filter, such as the 
median filter, the operation can NOT be performed in terms of convolution due the linearity of convolution). 
Fortunately, there exist some relationship between the convolution and Fourier transform as follows: 


S(g@%, W@f(% y)) = Glu, v)Flu,v), (3.7.1) 


where &! denotes the convolution operator ,.3 represents the Fourier transform operation. Here the Fourier 


transform of the result of convolution between two functions is the product of the Fourier transform of 
those two functions respectively. 


The Fourier transform provides a direct theoretical explanation of the filtering technique used in the 
spatial domain. Using this technique, we have an intuitive perception of which frequency will be used to 
filter the signal. A straightforward step in the filtering design is to select a filter which removes some 
signals at certain frequency bands that we do not want to keep. The design of the filter will be down to 
how to specify F (u, v) -- the Fourier transform of the spatial filter mask f(x, y). F(u, v) can be referred 
to as the filter transfer function. 
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A good example is the commonly used filtering mask -- Gaussian filtering. Figure 29 shows an Gaussian 
filter function in the space domain and its power spectrum in frequency domain. Note the difference 
between the variance. 


| 
1015 4 
| 

1 


(a) (b) 


Figure 29 A Gaussian filter in the frequency domain (a) and its power 
spectrum in the frequency domain (b). 
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3.8 Image Smoothing Using Frequency Domain Filters 


In this section we will introduce three main low pass filters including ideal lowpass Filters, Butterworth 
lowpass filters as well as Gaussian lowpass filters. In the following section, we will use the panda image 
to show the principles of the different filters: 


(a) (b) 


Figure 30 An image of a panda (a), and its power spectrum of the gray level image (b). 
3.8.1 Ideal Lowpass Filters 


An ideal low pass filter is to switch off some high frequency signals and to allow the low frequency signal 
to get through based on a threshold. Mathematically, it is defined as follows: 


0, u2+v2>T0 


8.1 
1, u2+v2 <T0 Ost) 


M(u,v) =| 


Here, TO is a threshold, usually a nonnegative number. It represents the square of the radius distance from the 
point (u, v) to the center of the filter, usually called cut-off frequency. Intuitively, the set of points at u? + 

v? = Tp forms a circle. It blocks signals outside the circle by multiplying 0 and passes the signals outside the 
circle by multiplying 1. By doing this, we get the following filtered results in the frequency domain: 


0, u2+v2>T, 


iu, v).*M (u,v) = i WD), wey < Ty 


(3.8.2) 


The operator .* represents the element wise dot product operation. Figure 31 shows an image example 
and its filtered image using the ideal low pass filters. 
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(a) (b) 


Figure 31 An ideal low pass filtering with cut-off frequency set as 10 (a), and the corresponding 
filtered image (b) from Figure 30 (a). 


3.8.2 Butterworth Lowpass Filter 
The ideal low pass filter has a sharp discontinuity as Ty. This will introduce some ring effects on the 


filtered image as shown in Figure 31 (b). In contrast, the Butterworth lowpass filter is much smoother at 
Ty and defined as follows: 


M(u.v) = ———> (3.8.3) 


T(u,v) _ 
1+| To | 
where T(u, v) is the radius from the point (u, v) to the origin in the frequency domain. 
Figure 32 shows an example of Butterworth low pass filter with Tg = 10, and the filtered image from 


Figure 29 (a). Comparing the filtered image with the one in Figure 30, we notice that the Butterworth low 
pass filter produces less ringing effects than an ideal low pass filter using the same cut-off frequency. 


(a) (b) 


Figure 32 Butterworth low pass with cut-off frequency set as 10 (a), and the corresponding 
filtered image (b) from Figure 31 (a). 
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3.8.3 Gaussian Lowpass Filters 


In 2D image processing, an exemplar filter is Gaussian filtering. It is perhaps the most commonly used 
low pass filter in the literature. Mathematically in the frequency domain, it is defined as follows: 


G(u,v) = Ae(—d(u, v)/207) (3.8.4) 


where d is an affine distance function from the center point to the current point (u, v), usually denoted by 
d= = (u? + v2) with o denoting the standard deviation. Figure 33 shows the plot of the 2D Gaussian 


low pass filter and the filtered panda image respectively. 


(a) (b) 


Figure 33 Gaussian low pass filter (a) and the smoothed panda image (b), with o set as 10 and 
Ais set to 1. 
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3.9 Image Sharpening Using Frequency Domain Filters 


3.9.1 Ideal Highpass Filters 


An ideal highpass filter is to switch off some low frequency signals and to allow the high frequency signal 
to pass based on a threshold. It implies the difference between 1 and the ideal low pass filtering. 
Mathematically, it is defined as follows: 


1, u2+v? =>T0 
M(u,v) = , - 3.9.1 
HY ty 402 <T0 ve 
In contrast to the ideal low pass filter, it blocks signals falling inside the circle by multiplying 0 and pass 
the signals settling outside the circle by multiplying 1. By doing this, we get the following filtered results 
in the frequency domain: 


F(u,v), u2+v2=Tp 


3.9.2 
0, u2+v2<Ty vee) 


F(u,v).* M(u,v) = 


Figure 34 shows an example of using the ideal highpass filters applied to an image. 


(a) (b) 
Figure 34 An ideal high pass filter (a) and corresponding sharpened image (b). 
3.9.2 Butterworth Highpass Filters 


The Butterworth Highpass Filter is simply the difference between | and the Butterworth lowpass filter T. 
It is defined as follows: 


M(u.v) = 1- (3.9.3) 


—— 
aco 


where again T(u, v) is the radius from the point (u, v) to the origin in the frequency domain. 
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(a) (b) 
Figure 35 Butterworth high pass filter (a) and the results of filtering the panda image (b). 
3.9.3 Gaussian High Pass Filters 
A Gaussian high pass filter is defined as follows: 
G(u,v) = 1— Ae(—d(u, v)/207), (3.9.4) 


: : : . 1 
where d is the distance from the center point to the current point (u, v), usually denoted by d = a (u? + 


v*) with 6 representing the standard deviation. 


Figure 36 shows a Gaussian high pass filter and the results of applying this filter onto the panda image. 


(a) (b) 


Figure 36 Gaussian high pass filter (a) and the results of filtering the panda image (b). 
3.9.4 The Fourier transform of Laplacian operator 


Another common high pass filer is the well known Laplacian filter. It is a 2D isotropic measure of the 
second spatial derivatives of an image. It can capture the region with rapid intensity changes. Traditionally 
it is used for edge detection. But in recent years it has been widely used for blob detection based on its 
extended version the Laplacian of Gaussian (LoG). It is mathematically defined as follows: 
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071 , O71 
V71(x, y) = ax2 + ay? (3.9.5) 


It is easy to prove that the Fourier transform of the above operator is: 

F(V7I (x, y)) = —4n(u2 + v2)I(u, v) (3.9.6) 
Thus the filter is as follows: 

F= —4n(u? + v?) (3.9.7) 
The Laplacian filter is sensitive to noise. In practical, it is common to blur the image first and then apply 
the Laplacian operator onto the smoothed image. The transform is traditionally termed as Laplacian of 


Gaussian (LoG). The Fourier transform of the LoG transform can be approximately represented as the 
following forms: 


—4(u2 + v?)exp (—0.5(u? + v?)/a7) (3.9.8) 
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(a) (b) 
' 
- . - 
al 
(c) 


(d) 
Figure 37 Laplacian fitler (up-left), the results of laplacian fitlered image (up-right). 


LoG filter (down-left). For better visibilty, here we show the inverse of the LoG filter, 
and the results of Log filter image (down-left). 


3.9.5 Unsharp filter, Highboost Filtering, and High-Frequency-Emphasis Filtering 
The unsharp filter is a very simple operator that subtracts a smoothed (unsharped) version of an image. 
This operation is called unsharp filter and the resulting difference is usually called unsharp mask. It is 
majorly used in the printing and publishing industry for enhancing image edges. 

In principle the unsharping mask is an edge image e(x, y) from the input image / (x, y) by the unsharping filter: 


ey) = 1G, y)— 1067) (3.9.9) 


where I(x, y) is a blurred version of the image I(x, y). In order to enhance or boost the edge image (high 
frequency part), the edge image e(x, y) is then added back to the original image a weighting strategy as follows: 


g(x,y) = I(x, y) + w * e(x,y) (3.9.10) 


Generally, w is non negative and it controls the contribution of the edge component of the image. It is can 
be easily seen that when w is set to be 0, there is no additional contribution of the edge part. When w is 
larger than 1, it is sometimes called high boost filtering in the spatial domain. 


We can easily perform the above operation in the frequency domain by applying the Fourier transform on 
the above equation. The final form can be as follows: 


G(u,v) = I(u,v) +w(1 — A(u, v)I(u, v)) (3.9.11) 
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where I(u, v) is the Fourier transform of the original image I(x, y) and G(u, v) is the Fourier transform of 
the boosted image g(x, y). 


3.9.6 Frequency Domain Homomorphic Filtering 


Image can also be expressed as the product of illumination and reflectance. This is known as the 
illumination-reflectance model. This model is first developed by MIT, USA. This model basically states 
that an image is produced by a sensor based on the amount of energy that it received from the object. The 
energy that an image sensor or camera received is determined by illumination radiation source and the 
features of the reflectance of the object itself. Their relationship is viewed as a product. It is 
mathematically expressed as follows: 


f@,y) =i, y) «r(x, y) (3.9.12) 
where f(x,y) is the image. i(x, y) is the illumination component and r(x, y) the reflectance component. 


The reflectance is the features of the object itself. Generally we are more interested in reflectance than 
illumination. A basic idea is to enhance reflectance and to reduce illumination. In order to do this, a common 
approach is called Homomorphic Filtering. To utilize the advantages of the Fourier transform, we need to 
convert the product into sum by applying the log operator on both sides of the equation. Thus we get: 


log(f) = log(i) + log (r) (3.9.13) 


Now the illumination-reflectance model in the log space becomes a linear model. It is easier for us to 
process this in the frequency domain. It has been shown that the illumination component majorly usually 
lies in the low frequency band and reflectance in the higher frequency band. This motivates us to borrow 
some high pass filter to enhance the reflectance part and degrade the other. Once doing this, we then can 
use the inverse Fourier transform to convert it back to the log space, and finally deploy the exponent 
operation to convert the log space into the original space. A common flow chart of the Homomorphic 
Filtering is shown in Figure 38. 


Figure 38 The flowchart of the homomorphic filtering procedure. 
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Filtering in the Frequency Domain 


Figure 39 shows an example of image enhancement using the homomorphic filtering. 


(a) 


(b) 


Figure 39 An example of image illumination correction using the homomorphic filtering. 
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Summary 


The Fourier transform is a key component of image processing. In other words, without Fourier transform, 
modern digital image processing will lack amounts of contributions to the community. In this chapter, we 
have introduced the principles of the discrete Fourier transform in 1D and 2D cases. We also illustrated 
the sampling theorem. Keep in mind that the sampling rate of a signal should be at least two times higher 
than the highest frequency of the image content. We further presented several low pass and high pass 
filters in the frequency domain including ideal low pass filter, Butterworth low pass filter, as well as 
Gaussian filter. A simple rule is that the sum of the low pass filter and the corresponding filter of the same 
type are usually equivalent to one. Other filters are also common in the field of image processing such as 
the median filter. However, those filters cannot be easily interpreted in the frequency domain. We have 
introduced them in Chapter 2troduce them in the next chapter. 
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Problems 


(12) What is Nyquist-Shannon sampling theorem? 
(13) W hat is the relationship of lowpass filter and high pass filter? 
(14) Explain the different steps and motivation of the homomorphic filtering. 
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4. Image Restoration 


4.1 Image degradation and restoration 


A degraded image is a clean image with an additive noise term. In a mathematical language, given a clean 
image f(x, y) and a noise function n(x, y), the degraded image is generated as follows: 


g(x,y) =A(x,y)* f(x,y) + n(x, y) (4.1.1) 


where A(x, y) is the transformation function from the clean image to the degraded image in spatial domain 
and “*” is convolution. 


Restoration is the converse process of degradation. That is to say, the purpose of image restoration is to 


seek a maximum similarity of the original clean image after an appropriate process is applied. Such a 
degradation/restoration procedure can be represented in the following graph: 


g(x,y) 


f(x,y) 


Degradation 
h(x, y) 


Restoration 


It is essential to understand how noise affects images before an optimal restoration method can be 
introduced. Without this knowledge, it will not be possible to discover a proper solution. In the nest 
section, the noise models will be introduced, where up to date research outcomes are presented. 


4.2 Noise analysis 


Image noise is normally due to environmental conditions (e.g. light, temperature, etc), sensor quality and 
human interference. Of these factors, the first two seems to be dominant in real practice. For example, a 
camera of low resolution easily leads to fuzzy imaging outcomes. A moving camera can collect image 
sequences with strong motion blurred effects. 


To perform appropriate noise analysis one of the possible ways is to understand the noise models behind 
the image signals. After these noise models have been identified, it may become much easier to apply 
corresponding strategies for problem solution. In general, these noise models can be grouped into time and 
frequency domains, the latter of which considers the difference between the frequency properties of the 
clean image and noise. More details can be found in the following. 


The first noise model is Gaussian. This is a well recognised noise model due to its mathematical tractability 
and popularity of use. To better understand the principle of Gaussian models, one of the options is to extract 
its probability distribution function (PDF). The PDF of Gaussian models is expressed by 


Download free ebooks at bookboon.com 


58 


Digital Image Processing — Part I Image Restoration 


Please click the advert 


a 
p(s)= ae il (4.2.1) 


where s is gray level, is the mean of s and a is its variance. The difference between the image signals 
and the average value is anti-proportional to their probability values. 


The second one is uniform noise model. The PDF of uniform noise is 


1 
p(s)=4a-b (4.2.2) 
0 


where the upper case results from a > s > b; otherwise the lower case stands. This indicates that the 
probability values will be the same no matter where each pixel locates. 


The third one is Rayleigh noise, whose PDF is delineated by 


_(s-ay 


2 b 
au ~A)JEAD (4.2.3) 
0 


P(s)= 
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where the upper case results from s > a (a and b are two constants); otherwise the lower one is valid. This 
noise model results a deformed shape of thr Gaussian model. 


Another quite common noise model is salt and pepper noise which has a PDF as follows 


P 


P(s)=4F, (4.2.4) 
0 


where the upper case is due to s = a, the middle case is due to s = b and otherwise the lower case is valid. 


Figure 40 illustrates some images with additive noise (the noise power is also included in the caption). 


Figure 40 Noisy images with different additive noise: (a) original image, (b) Gaussian (0, 0.2) 
(mean, variance), (c) speckle (0.2) and (d) salt & pepper (0.2). 
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4.3 Restoration with spatial analysis 


A number of approaches have been established to handle the image restoration problem. These normally 
consist of spatial and frequency based approaches. In this subsection only spatial analysis based methods 
are concerned. Of the developed schemes, spatial filtering is one of the frequently investigated areas and 
its applications have been well developed. Some examples are summarized here. 


Geometric mean filter: This is a model that incorporates geometric mean. Each restored image pixel is 
the mean value of the pixels falling in the investigated area. The underlying mathematical expression for 
the restoration is followed: 


fa@y=(T1f@om (43.1) 


(q,t)<O 


where (m, n) indicates the size of a subimage window. Examples of applying this filter to the noisy images 
(shown in Figure 40) are given in Figure 41. It demonstrates that this filter is more suitable to handle the 
case of speckle noise. 


Figure 41 Examples of using the geometric mean filter, where (a) Gaussian, (b) speckle and (c) 
salt & pepper. 


Median filter: This is a statistics based filtering algorithm. Any image pixel can be replaced by the 
median of the image pixels in the neighborhood of that pixel: 


f (x,y) 7 median f(q,t)} (4.3.2) 


Figure 42 illustrates the outcomes of using the median filter to the images shown in Figure 40. It is clear 
that this median filter has the best performance when we restore the salt & pepper image. 
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Figure 42 Outcomes of applying the median filter to the noisy images. 


Max and min filters: The median filter takes the mid point of the image intensity. On contrary, max and 
min filters are designed to extract the maximum or minimum point in the ranked numbers. This filter 
effectively works in some extreme cases, where the noise makes the image intensities outstanding. 
However, it will not perform well in the presence of regular noise distributions. 


Adaptive filters: These filters are intended to adapt their behaviors depending on the characteristics of the 
image area to be filtered. This is based on a simple fact that many signals actually are the combination of 
the known or anonymous noise resources and hence it is unlikely to obtain an appropriate noise model to 
describe the characteristics of the noise. 
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One commonly used adaptive filters is local noise reduction filter, where the denoising result can be 
represented by 


FO) =Fy)~ SU sy) — ml) (4.3.3) 


where o,” is estimated variance of the noise corrupting the image, m and o,’ are the mean and variance of 
the pixels in the sub-image. These variables sometimes can only be available through historic estimation. 
The performance of this filter is revealed in Figure 43. 


Figure 43 Illustration of the adaptive filter in image denoising. 
Wiener filter: This is a statistical approach to pursue an outcome that can minimize the mean square 


error between the restored and original images, given that the image and noise are uncorrelated with 
normal distribution: 


e = ELf(x,y)- f(x y)P (4.3.4) 


This leads to the outcomes shown in Figure 44. This shows that the restored images of speckle and salt & 
pepper noise are quite similar. 


Figure 44 Performance of the Wiener filter in image restoration. 
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4.4 Restoration with frequency analysis 
In the frequency domain the restored image can be represented as follows: 
G(u,v) = H(u,v)* F(u,v)+ N(u,v) 


Where G, H, F and N are the Fourier transforms of g, h, fand n, respectively. H is the optical 
transformation function. 


In the frequency domain the degradation process is equivalent to that of convolving the image with an 
optical transformation function. Similarly, the restoration can be referred to as deconvolution: 


N(u,v) 


F(u,v) = F(u,v)+ Hos) 


(4.4.2) 


This deconvolution process can be directly applied to solving the practical problems. As an example, the 
Wiener filter is particularly used here for further discussion. The Fourier transform of the restoration is 


F(u,v) 


_ H(u,v)S, (u,v) 
| S(u,v)| H(u,v)[? +55 (u,v) 


[eum 


| 1 | H(u,v)| 
| A(u,v) | H(u,v) [2 +S, (u,v)/S, (u,v) 


[eum 


where H(u,v) is degradation function, H’(u,v) is the complex conjugate of H(u,v), |H(u,v)|” = H' (u,v) 
H(u,v), S;(u,v) = |N(u,v)| that is the power spectrum of the noise, and S(u,v) = |F(u,v)|° that is the power 
spectrum of the clean image. The above equation indicates that the power spectra of the restored image 
depends on the power spectra of the transformation function and the ratio of the power spectrum of the 
noise and the clean image. 


The second example is notch filter. This filter has a similar concept to that of bandpass and Butterworth 
filters. It rejects or passes frequencies in a pre-defined areas around the central frequency. Now, the 
rejection case is investigated as an example. Let the centre and its symmetry be (wo, vo) and (-wo, -Vo) 
respectively. If the noth reject filter is applied to a radius D,, then the transfer function can be expressed as 


0 
H(u,v) = i (4.4.4) 
where D, = D; or D. < D2 leads to the upper case with 


D(uv) =[U-M/2-ujy+W-N/2-v/P 
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and 


D,(u,v) =[U-M/24+uy+WV—-N/24+v/P (4.4.6) 


The filter will look at an area with the shifted centre to (//2, N/2). Similarly, the transfer function of a 
Butterworth notch reject filter of order 2 is derived as follows: 


Ht oe (4.4.7) 
D; 2 
1+[ £ ] 
D,(u,v)D, (u,v) 
A Gaussian notch reject filter has the transfer function 
1 Pias)P2 (wv), 
Hiu,v)=1-e ” (4.4.8) 


Another typical example is constrained least squares filter. The Wiener filter has an evident drawback: 
It is desirable to know the power spectra of the clean image and noise. In spite of the existence of an 
approximation of the power spectra, S; and S, may not be constant and stable. Therefore, it is desirable to 
investigate the minmisation of the second derivative of the image for satisfaction of smoothness: 
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in| SS aes vr (4.4.9) 


x=0 y=0 
Subject to the following constraint: 
ll g— AF I?ll? (4.4.10) 


where ||~||? is the Euclidean vector norm, and V” is the Laplacian operator. The frequency domain 
solution to this minimization problem is given by the following equation: 


A * (u,v) 
| Hu,v) | +7 | P@,v) |? 


F(u,v) -| aan (4.4.11) 


where y is a parameter that can be tuned to satisfy the constraint equation. P(u,v) is the Fourier transform 
of the function: 


0 -l 0 
p(x,y)=|-1 4 -1 (4.4.12) 
0 -l O 


One of the optimisation techniques for computing y is to iterate the following procedure: 
Let the residual of the difference ( g — Hf ) be r;. We use a function of y to describe r, 

IN=rr Aly (4.4.13) 
Let 

II 7, IIa |? te (4.4.14) 


where c is an variable that controls the estimate accuracy. Once c is defined, an optimal y can be obtained 
by continuous adjusting in order to satisfy Equations 4.4.6 and 4.4.9. 


To compute |[r,||7, one can look at its Fourier transform: 
|R, P=Gw,v)-H(u,v)F (u,v) (4.4.15) 


Since 
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In t=¥ dr 
a= »V) (4.4.16) 


x=0 y=0 


Therefore, the variance of the noise over the entire image can be expressed as follows: 


M-\N-1 


a= Dnt y)-m,)° (4.4.17) 


x=0 y=0 


Where 


M-1\N-1 


m,, = ney) (4.4.18) 


x=0 y=0 
Eventually, the following equations stand: 
|| 2 ||’ = MN(o2 +m?) (4.4.19) 


This equation shows that image restoration can be performed if the variance and mean of the noise buried 
in the image are known beforehand. 


Figure 45 illustrates the performance of the contrainted least sqaures filter in image restoration, where the 
overall noisy images have been restored well. 


Figure 45 Examples of image restoration using the least squares filter. 


4.5 Motion blur and image restoration 


Motion blur is due to the relevant motion between the image sensor and the objects to be acquired during 
the shot period. Normally, this happens in the case where any motion is longer than the period of exposure 
constrained by the shutter speed. If this ocurres, the objects moving with respect to the sensor (i.e. video 
camera) will look blurred or smeared along the direction of the motion. 
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Motion blur can lead to special effects in computer graphics and visualisation. These effects sometimes 
are valuable. For example, in sports motion blur is used to illustrate the speed, where a slow Suter speed is 
normally taken. Figure 46 clearly shows motion blur in the areas of the running path and lower bodies. In 
the meantime, motion Blur has been Commonly used in the animation so as to create a synthetic scene for 
visualisation. One of the examples is shown in Figure 47. 


Figure 46 Motion blur in the racing scenario ([6)). 
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Figure 47 Motion blur in the train track scenario ([7]). 


However, motion blur also causes certain side-effects. One of these effects is the nuclear details in the 
blurred image areas. These areas need to be sharpened so that further process can be carried out, e.g. 3- 
dimensional construction that requires very accurate correspondences across the image frames. 


Before the blurred images can be deblurred, it will be very helpful to generate an optimal scheme if one 
understands how these images have been blurred. Figure 48 denotes a synthetic motion blurred image. 
The strategy used to create this image follows: 


e Move image pixels along a certain angle, /,; 
e Multiply the clean image / with the changed image /, in frequency domain, ,,; 
e Inverse Fourier transform of J,,. 


Figure 48 Illustration of a synthetic blurred image. 
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The methods used to deblur the blurred images can be those that have been introduced in the previous 
sections, e.g. mean, Wiener and constraint least squares filters. Figure 49 illustrates the deblurred image 
of Figure 48 using the Wiener filter. 


Figure 49 Illustration of the deblurred image using the Wiener filter. 
4.6 Geometric transformation 


The restoration techniques introduced in the previous sections are normally used to handle the images that 
noise signals are additive to. In this section, the case of changing the spatial relationship of the image 
points will be discussed. There are mainly two groups of geometric transformation: (1) Spatial 
transformation where the image pixels are re-arranged, and (2) image interpolation/extrapolation where a 
part of the images is filled or taken away. 


First of all, spatial transformation is introduced. Assume an image pixel has coordinates (x, vy). Then its 
gemetric distortion will have different coordinates as (x’, v’). The relationship between the changed 
coordinates is 


= Ly) (4.6.1) 


‘= T(x, y) 


S 
| 


Where 7; and 7> are two spatial transformation functions that generate the geometric changes. For 
example, a commonly used geomtric distortion is described in a pair of bilinear equations: 


me = Cx + Coy e Oxy + Cy (4 6 2) 


T,(%,y) = c.x + Coy + OXY + Cg 
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Secondly, image interpolation/extrapolation is discussed. Due to Eq. (4.6.2), the distorted image 
coordinates maybe non-integer. However, image pixels can only be integer. Therefore, some “pixels” will 
be interpolated or extrapolated. One of the common approaches is to use a nearest neighbor approach. 
This approach allows the inserted or extracted image patches to have the same intensity/colour as the 
closest pixels. For the simplicity purpose, Figure 50 shows a rotated image and its recovered image 
according to the estimated scale and angles. 


Figure 50 Illustration of the geometric distortion image (a) and its correction (b). 


Summary 


In this chapter, the concepts of image degradation and restoration have been introduced. The noise models 
are investigated before any image restoration algorithm is presented. These models consist of Gaussian, 
uniform, Rayleigh and salt & pepper noise. Afterwards, the corresponding techniques to restore the 
degraded images were also summarized. These techniques are catergorised into the spatial and frequency 
domains. The latter include the Wiener, notch and constrained least squares filters, while the former has 
examples such as the geometric mean, median, max and min filters, Wiener filter. Due to the importance 
and common use motion blur has been introduced in an independent section and deblur approaches are 
also summarized. Afterwards, geometric transformation is presented, where spatial transformation and 
image interpolation/extrapolation are individually described. 


In general, image restoration and the developments in this field are very important in the domain of image 
processing. The quality of the restored images will directly affect their applications and uses in the 
community. Some of the possible resources of image degradation and restoration have been presented. 
However, many new strategies are still in their infancy stages and hence more efforts will be made in 
order to move forward the developments. 
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Problems 


(15) 
(16) 
(17) 
(18) 
(19) 


(20) 


(21) 
(22) 
(23) 
(24) 


What is image degradation/restoration? 

Why do we need image restoration? 

Please summarize the noise models that have been presented in this section. 

How many groups of spatial analysis are available up to date? 

What is the principle of motion blur? Can you apply motion blurring to the image shown below? 


How can we handle the motion blurred images? Hints: Provide a general algorithm for the 
deblurring purpose. 

What is geometric transformation? How can you recover a spatial distortion? 

Can you use a different method for deblurring Figure 22? 

Please derive Equations (4.4.7) and (4.4.8). 

Please try to derive Equation (4.4.11). 
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